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* Department  of  Civil  Engineering ,  University  College 
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1.  INTRODUCTION 


Despite  a  continuing  effort  to  identify  optimal  finite  ele¬ 
ment  grids  most  of  the  finite  element  computations  today  still 
rely  on  an  a-priori  mesh  design  based  on  the  user's  intuition  and 
experience.  Once  the  mesh  is  designed,  however,  there  seems  to 
be  growing  evidence  that  high  order  isoparametric  elements  pro¬ 
vide  a  better  refinement  process  than  mesh  subdivision. 

- -^This  paper  is  concerned  with  the  identification  of  the  dis¬ 
cretization  error  in  finite  element  solution  and  the  definition 
of  optimal  refinement  processes.  The  advantages  and  limitations 
of  the  hierarchical  approach  f.;>8 ,  3,6 --4Q-]  will  be  discussed  and 
it  will  be  shown  how  the  intelligent  enrichment  of  the  finite 
element  grid  can  be  left  to  the  computer  if  a  capacity  for 
a-posteriori  error  estimation  exists  within  the  finite  element 
codey  [1-14,  37-40,  46-49] . 

2.  HIERARCHICAL  FINITE  ELEMENTS 

2. 1  Hierarchical  Shape  Functions 

The  concept  of  hierarchical  finite  elements  dates  from  1970 
and  these  were  first  introduced  [58]  with  the  objective  of  crea¬ 
ting  elements  that  would  allow  an  easy  transition  from  a  region 
where  a  finite  element  solution  required  a  high  degree  of  refine¬ 
ment  to  a  region  where  a  lower  degree  of  refinement  was  sufficient. 
Other  advantages  soon  become  apparent  and  it  will  be  shown  here 
that  the  hierarchical  concept  is  very  powerful  in  allowing  an 
error  indication  capability  that  can  be  used  for  adaptive  mesh 
refinement  [36-39,  46-49]. 

We  will  begin  by  defining  hierarchical  finite  elements  as 
those  in  which  successive  refinements  are  additive  in  the  manner 
of  additional  terms  in  a  Fourier  series.  It  follows  that  the 
"stiffness"  matrix  corresponding  to  the  hierarchical  element  at 
a  certain  level  of  refinement  is  a  sub-matrix  of  the  "stiffness" 


matrix  corresponding  to  a  higher  level  of  refinement. 

This  leads  to  matrix  approximation  equations  of  the  type: 
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where  (2.1)  is  the  finite  element  equilibrium  equation  correspon¬ 
ding  to  a  certain  formulation,  and  equation  (2.2)  is  the  same 
equation  corresponding  to  a  higher  order  or  refinement.  The 
matrices  and  q^  remain  unchanged. 

Consider  solving  the  linear  differential  equation 


A  (4>)  =  L  <!>  +  q  =  0  in  !1 

wi th  boundary  condition 

B  (cf>)  =  L  <{>  +  s  =  0  on  F 


(2.3) 


(2.4) 


where  T  is  boundary  of  f2 . 

^  Approximate  the  solution  $  of  (2.3)  (2.4)  by  the  solution 
4>  in  the  form 
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<J>  =  l  a  N 

,  mm 
m=l 


(2.3) 


with  a  proper  choice  of  the  basis  function  N.,  i  = 1 ,2, .  .  .  ,m. 

The  coefficient  a^  will  be  determined  from  the  condition 

that 

W  [L4>+q]dft  +  W  [L$+s]dr  =  0  (2.6) 

Q  6  _  T  e 

holds  for  all  Wg  and  Wfi,  where  We  and  We  are  functions  suit¬ 
able  for  the  problem  (2.3)  (2.4).  In  addition  (2.6)  is  meant  in 
a  generalized  way  so  that  it  can  be  used  even  if  LN^  does  not 
exist  in  a  classical  sense  (see  [59]). 

The  p-version  of  the  finite  element  method  is  defined  for 
a  sequence  of  solutions  as  the  enrichment  of  the  trial  and  test 
function  set  through  the  introduction  of  shape  functions  Nm  cor¬ 
responding  to  higher  order  polynomial  degree,  p,  while  the  h 
version  is  defined  as  the  approach  equivalent  to  the  reduction 
of  the  finite  element  mesh  size  h,  maintaining  constant  the 
polynomial  order  [8-9,  46-49]. 

In  1-D  the  p-version  of  the  finite  element  method  involves 
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the  addition 
functions  as 
expressed  as 
-1  <  £  <  1 
chical  shape 


to  the  linear  element  of  quadratic,  cubic,  etc.  trial 
shown  in  Figure  1.  Considering  the  trial  functions 
a  function  of  the  local  element  coordinates  with 
we  may  have  for  the  p-version  the  following  hierar- 
functions.  For  the  linear  terms 
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For  the  quadratic  the  'obvious'  form  will  be  a  quadratic  that 
goes  to  zero  at  points  £  =  -1  and  £  =  1  so  that  it  does  not  inter¬ 
fere  with  the  a  and  a,  coefficients.  Thus 
o  1 

N2  =(£  -1)  (£  +1)  (2.8) 


For  the  cubic  and  higher  order  elements  the  only  restriction  is 
that  the  shape  functions  will  go  to  zero  at  £  =  -1,  £  =  1  ,  so 
that  we  write  for  the  p-th  order 

Np  =  (£  -1)  (£  +1)  £P~2  ,  p  >  2  (2.9) 

The  above  set  is  obviously  not  unique  and  many  alternatives  are 
possible. 

A  convenient  form  of  hierarchical  functions  is  given  in 
Peano  et  al.  [38]  as 
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where  p  5  2  is  the  polynomial  order. 

It  is  easy  to  observe  now  the  associated  variables  have  the 
meaning  of  higher  derivatives  of  4>  i.e.. 
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Although  the  successive  importance  of  these  hierarchical 
variables  diminishes,  the  optimal  form  is  one  which  gives  an 
orthogonal  set  of  shape  functions  in  relation  to  the  energy  inner 
product.  Such  a  set  of  shape  functions  will  be  given  by  integrals 
of  the  Legendre  polynomials  in  the  following  form 
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if  the  differential  equation  -<j>  +  q  =  0  is  considered. 

This  set  was  introduced  by  Zienkiewicz  et  al.  [58].  Again, 
the  multipliers  have  the  meaning  of  a  measure  of  the  higher  deri¬ 
vatives  of  0  at  the  centre  of  the  finite  element.  The  advantage 
of  this  set  of  hierarchical  functions  is  that  for  1-D  problems 
the  coupling  between  different  higher  order  degrees  of  freedom  is 
non-existent.  We  will  discuss  this  point  further  in  Section  2.2. 
We  shall  also  see  that  this  set  of  hierarchical  functions  will 
play  a  very  important  role  in  the  error  analysis  study. 

For  the  h-version  there  is  also  an  infinite  number  of  hier¬ 
archical  refinement  possibilities  corresponding  to  the  sub-divi¬ 
sion  of  the  initial  element  in  equal  or  unequal  parts.  If  we 
consider  a  hierarchical  refinement  corresponding  to  Figure  2,  we 
will  have  as  shape  functions 


N 


0 


1  -g 
2 


N 


1 


Lii 

2 


£  +1 

if 

£  <  o 

-e+i 

if 

£  >  o 

(2.13) 


etc.  Considering  Figure  2  it  is  obvious  that  the  physical  meaning 
of  the  linear  hierarchical  variables  is  a  relative  displacement 
set . 

Once  the  one-dimensional  interpolation  formulae  have  been 
established  the  generation  of  hierarchical  shape  functions  for 
rectangular  elements  is  almost  trivial  as, 

a)  the  corner  node  functions  are  simply  bi-linear  products,  and 

b)  'hierarchical'  functions  of  the  type  defined  above  are  always 
zero  at  the  corner  nodes. 


Polynomial  shape  functions  of  all  orders  in  two  dimensions  can  be 
obtained  by  simple  products  of  the  one-dimensional  formulae,  but 
in  general  losing  the  properties  of  the  orthogonality  mentioned 
above.  The  identity  of  hierarchical  variables  on  any  element  side 
with  those  on  the  adjacent  element  then  automatically  guarantees 
the  uniqueness  of  the  polynomial  along  that  side. 

The  three-dimensional  case  and  the  triangular  based  finite 
elements  are  just  special  cases  of  the  concepts  expanded  above 
and  we  direct  the  interested  reader  to  [36,  62]. 

The  hierarchical  and  non-hierarchical  shape  functions  for 
the  h  and  p  versions  of  the  finite  element  method  are  presented 
in  Figure  3  and  Figure  d  for  two-dimensional  problems.  It  is  easy 
to  show  that  a  direct  transformation  "rom  the  hierarchical  to  the 
non-hierarchical  formulation  and  vice-versa  is  possible.  For 
example,  it  is  possible  to  transform  a  quadratic  serendipity 
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element  to  hierarchical  form  by  a  matrix  transformation  after 
the  non-hierarchical  matrix  has  been  determined.  It  should  also 
be  noted  that  when  the  refinement  is  to  the  same  order  in  Figures 
1  to  4,  the  finite  element  solution  <p  in  (2.5)  from  the  hierar¬ 
chic  and  non-hierarchic  formulations  will  be  identical. 


2.2  Advantages  of  Hierarchical  Formulations 


In  1-D  the  advantages  of  a  hierarchical  formulation  over  a 
non-hierarchical  one  are  obvious  because  we  can  obtain  a  system 
of  equations  of  the  form 
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where  are  diagonal  matrices  if  the  functions  are  of  a 

suitable  orthogonal  form.  This  implies  an  improved  conditioning 
of  the  assembled  system  of  equilibrium  equations,  compared  to  the 
non-hierarchical  formulation  which  does  not  have  the  same  strong 
diagonal  character  and,  secondly,  the  possibility  of  a  direct 
solution  for  the  hierarchical  variables. 

In  2-D  we  expect  the  same  advantages  to  hold  although  the 
system  of  equations  is  now  not  completely  orthogonal  and  has  the 
form 


rKn 


L.K21 


(2.15) 


The  first  advantage  noted  above  is  carried  to  the  two- 
dimensional  case  by  the  fact  that  for  each  element  the  condition¬ 
ing  of  the  hierarchical  stiffness  matrix  is  better  than  the 
conditioning  of  the  non-hierarchical  one.  The  off-diagonal  links 
are  weakened,  implying  an  overall  better  conditioning  of  the 
resulting  system  of  equations.  This  is  in  fact  observed  by  num¬ 
erous  authors  [30,  53-56],  who  have  proved  that  a  relative  dis¬ 
placement  formulation  of  the  finite  element  method  yields  better 
conditioned  matrices  than  the  classic  total  variable  approach. 
Indeed,  the  same  can  be  said  of  so-called  local-global  element 
f  orms . 
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In  relation  to  the  second  advantage  we  can  say  that  the 
hierarchical  formulation  is  optimal  because  it  allows  for  all  the 
information  to  be  passed  from  one  discretization  level  to  the 
second  discretization  level  once  a  mesh  refinement  is  decided. 
Also  the  implicit  substructuring  existent  in  this  multi-level 
formulation  allows  for  the  very  effective  use  of  block  iteration 
solution  schemes.  This,  associated  with  the  better  conditioning 
ot  the  overall  system,  will  imply  a  fast  rate  of  convergence  for 
the  iterative  equation  solver.  Wachspress  two-level  elements 
[54]  are  another  way  of  achieving  these  objectives  although  on  a 
two-level  theory,  rather  than  a  multi-level  one.  See  also  the 
local-global  formulation  of  Mote  [ 30 ]  as  another  possible  form 
of  hierarchical  formulation. 


3.  ERROR  ANALYSIS  (A  HIERARCHICAL  APPROACH) 

A 

The  function  <{>  of  the  form  (2.5)  is  not  in  general  the 
exact  solution  of  the  problem.  Therefore  we  cannot,  in  general, 
have 

A 

L  4>  +  q  =  r  =  0  (3.1) 

L  4>  +  s  =  p  E  0  (3.2) 

(if  (3.1)  (3.2)  would  be  satisfied,  we  would  have  0  =<J>). 

Assume  for  simplicity  that  (3.2)  (i.e.  0=0)  holds;  then 
r  f  0.  Function  r  usually  (e.g.  if  L  is  a  second-order  equa¬ 
tion  and  have  not  continuous  derivative)  can  be  written  in 

the  form 


r  =  r][  +  r2  (3.3) 

where  r^  is  the  (usual)  regular  part  of  the  residuum  inside 
every  element  and  r0  is  the  singular  part  (Dirac)  function  con¬ 
centrated  on  the  interface  between  elements  with  physical  inter¬ 
pretation  of  concentrated  forces  whose  origin  is  indicated  in 
Figure  5. 

Considering  e.g.  equation 

-<t>"  +  q  =  0  (0)  4»(L)  =0  (3.4) 

and  using  piecewise  linear  elements  (as  in  Figure  5)  we  get  the 
regular  part 

rl  = 

and  the  singular  part  r^  are  the  forces  (Dirac  function)  concen¬ 
trated  at  the  nodal  point  with  the  magnitude  of  the  jump  of  the 
derivative  of  the  approximate  solution  on  the  given  nodal  point. 

For  the  potential  problem  in  two  dimensions  the  Dirac  function 
is  concentrated  on  the  element  interface  and  has  the  magnitude 
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J  =  .94.  - 


3n  ^  3n  2 


where  ^  .  i=l,2  is  the  normal  derivative  on  the  left  and  right 

side  of  the  interface.  In  the  case  of  elasticity,  Js  is  the 
traction  discontinuity  between  elements. 

3.1  Error  Definition  in  the  Energy  Norr. 

Denoting  e  =  <J>  -  <f> ,  the  error  of  the  finite  element  solution, 
then  obviously  subtracting  (2.3)  (3.1)  and  (2.4)  (3.2)  we  get 


Le  =  -  r 


L  e  =  -  p 


(3.5) 


and  we  assume  for  simplicity  as  before  that  P  =0.  The  goal  is 
now  to  measure  the  magnitude  of  e .  If  the  set  of  the  trial 
functions  and  test  functions  is  the  same,  then  for  a  suitable 
class  of  linear  problems  we  can  define 


=  |eLedO  -  =  £  j  C  -P)  L«-P  )d9.  '  (3.6) 


the  so-called  energy  norm,  and  through  it  we  can  measure  the 
magnitude  of  e . 


FIG. 5.  Interface  residuals  for  error  analysis. 
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where  II  *||„,T  .  is  the  norm  of  the  error  of  the  element  I-. 
"  11  E  ( 1  j  )  J 

This  process  obviously  needs  an  estimate  of  a£+^- 

In  one  dimension,  when  we  choose  as  shape  functions  poly¬ 
nomials  whose  derivatives  are  orthogonal,  we  have  from  (2.14) 
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qi  +  l 


K. 


(3.16) 


i+1 , i+1 


Here  q.  ,  can  be  determined  from 
l+l 


li+l 


=  N. 


i+1 


q  dil 


N. 


i+1 


(r-  L<t>)  dO  =  j  N^  +  ^  r  dO(3.17) 


because 
above . 


Ni  +  i  LO  dil  =0  due  to  the  orthogonality  referred  to 
Now  (3.11)  can  be  written  as 


Hell 


2  .  <W 

E<Ii>  EL+1 ,i+l 


(3.18) 


which  is  the  error  indicator  presented  in  [37,40]. 

When  the  shape  functions  do  not  possess  a  complete  orthogon¬ 
ality  a  degree  of  approximation  has  to  be  introduced  here. 

First,  has  to  be  estimated  using  the  previously  found  aj 

values  as 
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(3.20) 


because 


N.  . 

i  +  l 


dil  f  0 


(3.21) 


and  q^  +  ^  is  again  evaluated  from  the  previously  determined  . 

These  approximations  are  generally  tenable  if  near  orthogonality 
of  hierarchical  functions  exists, as  is  often  the  case  with  the 
p-type  elements. 

A  more  serious  shortcoming  of  this  error  estimate  is  how¬ 
ever  immediately  apparent.  If  the  residual  r  is  orthogonal  to 

N.  ,  in  (3.15)  the  error  indication  will  be  zero,  since 

i  +  l  ’ 
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j  N.+1-  r  dfl  =  0  (3.22) 

This  implies  that  the  proposed  criterion  leads  to  an  indication 
of  the  error  absorbed  by  a  hierarchical  refinement  on  the  exist¬ 
ing  mesh  rather  than  an  accurate  estimate  of  the  error  in  the 
finite  element  solution.  Therefore  we  will  denote  this  simply 
as  an  error  indicator  and  search  further  for  a  true  error 
"estimate" . 

3.3  Error  Estimation  in  One  Dimension 


Such  an  estimate  can  be  obtained  from  a  simple  calculation, 
when  e  =  4>  —  <j)  is  the  exact  response  to  r. 

Consider  the  model  problem 


+  q 


=  0 


(3.23) 


with  boundary  condition 


(f)(0)  =  4>(L)  =  0 

with  e  being  the  error  involved  in  a  linear  approximation  with 
nodal  points  x.  .  As  in  (3.5),  we  have 

d^e 

-  — j  +  r  =  0  with  e(0)  =  e(L)  =  0  (3.24) 

dx 


In  addition  it  can  be  shown  that  e(xj_)=0  in  our  particular  case 
(e.g.  [59]),  so  e  can  be  determined  on  every  interval  Ij  = 

( x i — i  j*;)  separately. 

J  J  —  77X 

Assuming  that  r=  rsin— ,  then  the  differential  equation 


d2e 


dx 


=  0 


e (0)  =  e (L)  =  0 


(3.25) 


is  easily  solvable. 


Substituting 


—  .  TTX 

e  =  e  sin— 


in  (3.25)  we  get 


TT  —  .  TTX  —  .  TTX 

~2  e  sm  —  =  -  r  sinr 

.  .  L2 

giving  e  = - 2  r 

7T  ** 

Further  on  L 


(3.26) 

(3.27) 
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So  far  we  assumed 
dual  will  involve  many 
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cnat  r  =  r  sm  -r-  . 
terms  and 


In  general,  the  resi- 
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It  can  be  easily  shown  that  we  get  (see  Appendix) 


II 


2  , 
r  dx 


(3.29) 


i.e.,  we  replace  equality  in  (3.28)  with  inequality  in  (3.29). 
Applying  now  (3.28)  to  every  interval  separately  [because 
e(x£)=0,  it  is  possible]  then  we  get 


vr 


where  *i+l-x.  =  h.+r 


(3.30) 


The  bounding  inequality  in  (3.30)  is  valid  for  r  of  all 
orders  of  polynomial  variation  in  the  one-dimensional  problem 
(3.23)  and  linear  elements.  Note,  however,  that  in  the  limit  of 
h  refinement  of  the  linear  elements  we  expect  that  r  will  be 
nearly  constant  on  each  element.  In  this  limit  only,  the  addition 
of  a  hierarchic  quadratic  term  will  give  a  valid  error  estimate. 
Here  we  can  rewrite  (3.14)  on  every  single  element  I  of  length  h 


<p  =  p  +  x(h-  x)  +  1 

and 

e  =  a.+1  x(h  -x) 

From  (3.24) 

r  =  -  2a .  , 
i  +  l 

and  substituting  (3.33)  in  (3.32) 


e 


-  j  x(h  -x) 


Then 


/ 

e  r  dl 


h 

( 

I  x(h-x) 
0 


dx  = 


r^h 

12 


3 


(3.31) 

(3.32) 

(3.33) 

(3.34) 

(3. 35) 
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Note  that  (3.32)  will  give  the  same  measure  of  the  error  for  con¬ 
stant  r  if  the  factor  is  changed  to  12.  We  therefore  imple¬ 

ment 

IMI l  =  iy  \  1+1r2  dx  (3.36) 

x. 

i 

as  the  asymptotically  correct  error  estimator  for  linear  elements 
and  one-dimensional  problems.  This  is  indeed  the  form  presented 
by  Babuska  et  al.  [4-6]  where  a  different  and  detailed  mathemati¬ 
cal  justification  is  given. 

3.4  Error  Indication  and  Estimation  in  Tuo  Dimensions 


In  two  dimensions  the  hierarchical  functions  have  their 
support  on  either  one  or  two  elements,  as  shown  in  Figure  6. 
Again  we  can  use  the  hierarchical  modes  to  make  a  total  projec¬ 
tion  as  in  (3.15).  We  note,  however,  that  in  two  dimensions  the 
hierarchical  contributions  to  K  are  not  diagonal  as  in  the  1-D 
equation,  and  the  approximation  of  (3.19)  has  to  be  used,  i.e. 


i+1 


K. 


(3.37) 


i+1 , i+1 


(for  j#  i+1  and  summation  convention)  then 

2 


II  e 1 1 . 


rq.  ,  +  K.  ,  .  a . 1 
lMi  +  l  l  +  l,  j 

K. 


(3.38) 


;i+l  i+1, i+1 

where  the  subscript  i+1  refers  to  the  new  hierarchic  mode. 


FIG. 6.  Hierarchical  supports  for  two-dimensional  error  analysis. 
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Note  that  if  we  sample  each  hierarchical  mode  independently 
(set  other  hierarchical  amplitudes  a^ =  0)  then  the  interaction 
of  the  new  modes  in  the  solution  is  ignored,  introducing  a  new 
approximation  in  addition  to  the  problem  of  this  error  indicator 
being  a  projection  of  the  true  error  in  the  next  hierarchical 
mode  and  not  the  total  error.  The  effect  of  this  new  approxima¬ 
tion  cannot  be  evaluated  a-priori. 

To  estimate  the  total  error  we  will  take  some  guidance  from 
the  hierarchic  nodes  in  Figure  6  but  attempt  to  base  the  analysis 
on  the  theory  developed  for  the  one-dimensional  case.  The  hier¬ 
archic  modes  indicate  that  not  only  must  we  consider  the  contri¬ 
bution  from  the  residual  on  the  element  (regular  part  of  the 
residual)  but  also  the  contribution  from  the  interface  (singular 
part  of  the  residual).  In  one  dimension  only  the  former  was  con¬ 
sidered  because  the  hierarchic  nodes  were  all  internal  to  the 
element . 

Consider  first  the  regular  part  of  the  residual  on  the  ele¬ 
ment.  We  will  over-estimate  the  error  associated  with  the  hier¬ 
archic  node  in  Figure  6a  by  releasing  boundary  conditions  to  give 
two  one-dimensional  responses  (see  Figure  7a)  and  associating 
half  the  residual  in  each  direction.  In  each  direction  (3.30) 


gives 


(3.39) 


so  that  the  addition  of  both  one-dimensional  contributions  gives 


r2dft  (3.40) 

‘J 

The  influence  of  the  singular  part  of  the  residual  can  be 
treated  in  the  same  way  if  the  residual  is  distributed  as  indi¬ 
cated  in  Figure  7b.  Now  we  have 


r  - 


b 

~h 


and  (3.30)  gives 


2  ,  - 

r  d  Q 


<Js 


l  h 

'  J2  d  y 
J  S 


dil 


0.41) 


O  .  42) 
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FIG. 7.  Deformation  patterns  for  two-dimensional  error  analysis. 


The  error  estimator  will  consider  each  element  separately. 

We  introduce  a  factor  of  1/ 2  into  (3.42)  because  each  boundary 
will  be  considered  twice.  The  error  estimate  becomes  for  the  bi¬ 
linear  element 


dfi 


■)  +  2 

2TT  \r 


O 

j“  dr 


(3.43) 


Since  this  derivation  is  based  on  a  superposition  of  one¬ 
dimensional  analyses  we  should  obtain  the  one-dimensional  estimate 
if  the  boundary  conditions  impose  essentially  a  one-dimensional 
problem  on  the  two-dimensional  mesh.  Using  (3.43)  and  constant 
right-hand  side  q  in  (2.3)  we  get  an  estimate  which  is  not  cora- 
plately  identical  with  the  one-dimensional  estimate.  Therefore 
we  reduce  the  estimate  by  modifying  (3.43)  to 


"  -"ec:.) 

j 


-  h  f  (  -)2  r 

1  — 7  !  I  r  -  r  J  d 

',T7~  J1.!. 

J 


h  f  2 
+  - y  I  _  J  d . 


(3.44) 


?T.  “  J 
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where  r  is  the  mean  value  of  the  residual  on  the  element.  This 
form  is  mathematically  justified  in  [2], 

It  follows  that  on  the  entire  domain  we  get 


in 

i 


(r  -  r  j  ^  dft 


+ 


(3.45) 


where  the  sum  is  taken  over  all  elements  and  I\  (if  Js  at 

the  boundary  is  properly  defined). 

We  derived  these  estimates  under  various  assumptions.  The 
question  arises  whether  they  are  acceptable.  The  answer  is  posi 
tive  (see  for  example  Babuska  et  al.  [2,7]  ).  There  the  term  1/  2 
is  replaced  by  ^  / 12  Co  obtain  the  asymptotically  correct 
estimate 


n. (r_r)2  dfi  + 


1 


(3.46) 


The  estimate  has  essentially  two  parts.  The  one  related  to 
the  regular  part  of  residual  (volume  integral)  and  the  other  one 
to  the  singular  part  of  the  residual  (jump  of  derivatives  term). 
It  can  be  theoretically  shown  that  the  first  term  (volume  inte¬ 
gral)  is  in  the  limit  negligible  with  respect  to  the  second  one. 
Practical  experience  shows  that  the  first  term  is  relatively 
small  also  for  coarse  meshes. 

The  estimate  is  called  a  svmpt  ot  ica  llv  correct 
when  the  ratio  of  the  right  and  left-hand  side  of  (3.46)  goes  to 
one  as  the  error  goes  to  zero.  The  estimate  (3.46)  has  this  pro¬ 
perty  when  some  mathematical  assumptions  are  satisfied.  One 
major  one  is  that  the  element  error  estimators  are  about  the  same 
ini  magnitude.  This  can  be  achieved  e.g.  by  an  adaptive  selection 
of  the  elements.  The  experience  shows  e.g.  that  the  asymptotical 
correctness  is  not  achieved  when  the  solution  has  singular  behav¬ 
iour  and  a  uniform  (obviously  improper)  mesh  is  used. 


3.6  Practical  Error  Analysis 

Both  the  estimators  (3.36)  and  (3.43)  and  the  indicators 
(3.18)  and  (3.38)  can  play  a  fundamental  role  in  the  finite  ele¬ 
ment  analysis.  The  estimators  allow  for  an  evaluation  of  the 
total  error  in  energy  in  the  current  finite  element  solution  and 
the  indicators  allow  a  rational  increase  of  degree  of  the  element 
(made  in  a  hierarchical  way)  and/or  element  subdivision. 

The  latter  is  a  direct  c  cnsequence  of  the  hierarchical 
error  indicators  being  a  projection  of  the  error  in  the  new  hier¬ 
archical  modes,  reflecting  thus  the  capacity  for  the  new  modes  to 
'absorb'  error  in  energy.  In  addition,  the  possibility  of  obtain¬ 
ing  an  accurate  estimate  of  the  error  in  an  appropriate  norm 
allows  the  program  to  stop  automatically  when  a  certain  accuracy 
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r 


has  been  achieved  or  to  indicate  the  order  of  accuracy  when  a 
certain  pre-specif ied  solution  cost  has  been  attained. 

It  seems,  however,  that  to  obtain  both  advantages  we  have 
not  only  to  compute  all  the  residuals  and  stress  discontinuities 
at  every  stage  of  the  interation  process  to  evaluate  the  esti¬ 
mate  (3.46),  but  also  all  the  hierarchical  stiffness  coefficients 
corresponding  to  the  possible  new  refinements  to  evaluate  the 
indicator  (3.38).  This  is  not  necessarily  the  case,  since  the 
error  indicators  can  be  obtained  as  projections  of  the  computed 
residuals  and  stress  jumps  in  the  new  hierarchical  modes,  i.e. 
in  (3.38) 


q .  +  K.  .  a. 
3  Ji  1 


r.  N.  dfl 
1  J 


+ 


N. 

J 


d$l  - 


N.  dr 
J 


where  anc*  T  i/2  and  Js,  ate  indicated  in  Figure  7b.  This 

computation  can  be  achieved  locally. 


3.5.1  Requirements  for  Practical  Error  Estimates .  We  consider 
that  a  practical  error  estimator  should  satisfy  the  following 
conditions : 


b. 


Be  determined  a-posteriori  from  information  defined  on  a 
local  basis .  , ,  ,  , 


If  we  define  an  effectivity  index 

the  ratio  of  the  predicted  energy 
exact  value  of  this  norm) ,  then  we 
meshes,  and  to  provide  reasonable 


0  = 


TO. 


(i.e. 


'exact 

norm  of  the  error  to  the 
require  9  >  1  for  all 
bounds  1  $  6  j-  2 . 


c.  Asymptotic  convergence  6 


d.  A  direct  interpretation  of  errors  in  stresses  should  be 
available . 


In  the  examples  that  follow  we  show  that  error  estimates 
which  satisfy  these  criteria  are  available  from  the  estimators 
of  the  form  (3.36)  and  (3.46)  if  the  foil  owing  amendments  are 
incorporated : 


1.  Theorem  2  of  [7]  states  that  there  exist  constants  k0  >  >  0 

independent  of  the  mesh,  such  that 


k,  II  ell, 


<  k, 


exac  t 


with  !!  ej|E  given  by  (3.4h).  Experience  has  shown  k-j  ;  0.3  and 
k2  1.6  on  a  large  number  of  problems,  and  both  are  asymptoti- 
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cally  equal  to  1  for  uniform  meshes  and  smooth  solutions.  To 
prevent  gross  violation  of  the  second  condition  required  of  the 
error  estimates  we  seek  a  value  of  k2  by  defining  a  factor  k2i 
for  element  i  as 


where 


k2i 


1  +  a 


f  u 

V**h 


(3.47) 


(3.48) 


is  the  energy  on  the  i-th  finite  element.  The  corrected  error 
estimate  becomes 


(3.49) 


We  have  taken  a  =2  in  all  applications. 

2.  For  isoparametric  transformation  of  the  element  there  is 
more  than  one  possible  choice  for  Che  length  parameter  h. 
Following  [19]  and  because  we  seek  an  over-estimate  of  the  error 
in  the  solution,  we  choose  h  as  the  length  of  the  maximum  side 
of  the  element. 

3.  Error  estimators  have  not  been  developed  for  elements  of 
higher  order  than  linear.  However,  significant  work  is  being 
done  in  this  area  by  Szabo  and  his  co-workers  [46-49]. 

The  sine  function  analysis  affords  the  following  extension. 

Quadratic  elements  would  match  the  predominant  part  of  the 
error  in  the  first  sine  function  i=  1,  used  for  linear  elements. 
Assuming  now  the  error  in  the  form  of  the  sine  function  with  i=2 
we  would  replace  the  coefficient 


(3.50) 


with  n=  8.  For  the  cubic  element  we  could  progress  to  the  mode 
i  =3  to  get  the  coefficient  with  n=  18. 

Thus  with  changing  polynomial  order  p  on  the  elements  the 
coefficient  of  the  estimator  can  be  replaced  by 


Tt  is  to  be  noted  however,  that  this  assumes  that  quadratic 
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interpolants ,  for  example,  can  completely  eliminate  error  in  the 
form  of  the  first  sine  function.  This  of  course  is  not  the  case 
so  the  factor  above  would  be  optimistic. 

Again,  in  the  interest  of  producing  an  over-estimate  of  the 
error  in  the  solution,  we  reduce  the  power  on  p  to  one  so  that 
a  general  form  of  (3.46)  for  elements  of  polynomial  order  p 
becomes 


dQ  + 


h 

24p 


Jr. 

1 


j2  dr 

s 


(3.51) 


Note  we  relax  the  requirement  9 ->■  1  here  in  2.  above. 


4.  The  analysis  above  strictly  applies  only  to  the  Laplace 
operator.  However,  it  can  be  generalized,  for  example,  to  pro¬ 
blems  governed  by  the  Navier  equations  of  elasticity.  Formula 
(3.51)  has  been  found  adequate  when  the  influence  of  Poisson's 
ratio  has  been  incorporated.  We  take 


!  I  e  1 1 


2 

E 


(n.) 


=  c 


an  + 


24p 


j2  dr  1 

r.  s  > 


with 


C  = 


1  -v 
E 


for  plane  stress. 


(3.52) 


5.  Local  estimates  of  the  error  in  stress  are  the  elusive  goal 
of  most  practical  error  analysis.  It  is  interesting  that  the 
error  estimates  advocated  here  are  evaluated  locally  element-wise 
are  justified  locally  by  appealing  to  St.Venant's  Principle  and 
are  backed  up  by  a  corrective  factor  which,  in  a  limited  number 
of  experimental  problems,  ensured  local  bounds  on  the  energy  of 
each  element. 

We  can  suggest  the  following:  a  bound  on  some  average  of 
the  stress  on  the  element  will  be  obtained  bv  scaling  local 
stresses  by  a  factor  FS . 

l 


FS. 

r 


(3.53) 


where  ||u  ||  is  the  prediction  of  the  energv  of  the  exact 

Eav 

solution  on  the  i-th  element  evaluated  as 


Ero 


(‘V 


' !  uh ! ' 2 

Lc..) 


r 


!  ^  i  ■  p  -  \  I 

E(  i’l 


(3.54) 
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and  1 1 u  1 1 

E(%) 


is  given  by  (3.48) . 


3.6  Applications  of  the  Error  Estimators  of  Section  3.5. 

3.6.1  Cantilever  with,  applied  end  moment,  linear  elements.  The 

simple  cantilever  shown  in  Figure  8  was  analysed  for  plane  stress 
using  the  four-node,  bilinear  element.  The  exact  solution  has  a 
quadratic  variation  of  vertical  displacement  along  the  beam. 

The  finite  elements  used  are  compatible  and  fully  integrated 
so  the  finite  element  solution  u^  is  known  to  be  stiff.  The 
error  estimate  ||e||g  does,  however,  over-estimate  the  error  in 
the  solution.  A  comparison  of  the  results  also  supports  the  con¬ 
stant  C  in  (3.52)  used  to  include  the  influence  of  Poisson's 
ratio  V  .  Finally,  the  moment  is  constant  along  the  beam  in  this 
simple  example  so  that  the  stresses  in  each  element  are  identical 
and  ||e||  from  (3.52)  is  the  same  for  each  element.  The 

E  (fl- ) 

error  estimates  in  the  table  can  therefore  be  interpreted  locally 
as  well  as  globally,  and  implementation  of  (3.53)  will  obviously 
give  good  estimates  of  the  local  error  in  stress. 

The  last  two  columns  in  the  table  give  the  ratios 


respectively.  This  measure  may  provide  the  best  qualitative  indi¬ 
cation  of  the  accuracy  of  the  finite  element  stresses. 

Note  that  all  calculations  between  values  in  these  tables 
must  be  based  on  the  formula 


exact 


This  formula  expresses  the  fact  that  the  error  is  orthogonal  to 
the  finite  element  solution. 


3.6.2.  Small  circular  hole,  h  and  p  refinement.  The  configura¬ 
tion  of  the  problem  and  the  finite  element  meshes  used  are  shown 
in  Figure  9.  The  analysis  was  for  plane  stress  with  v  =  0.3.  The 
results  given  as  Case  1  in  Table  2  are  for  the  bilinear  four-node 
element  and  subdivision  of  all  elements  into  four,  as  show  in 
the  figure.  In  this  analysis  the  surface  geometry  of  the  hole 
was  updated  in  the  refinement.  Again,  the  global  error  estimates 

!! e  !  I ^  are  exce 1  lent . 
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Initial  mesh 


TABLE  2 


Car  ae 

Plate  ■■■■nth 

Snail  Circular  Hole  . 

=  C.  I 

Case  1.  Linear  Elements 

(Complete  Subdivision) 

No.  Elements 

li  uh  'iE 

Error  Estimate 
i'  e  ;,£ 

0 

Error  Estimate 

!’ e  i;e 

-  * 

2- 

9t> 

"Exact" 

1.13981c 

1.  139829 

.583  x  10"- 
.393  x  10-2 

.85 

.81 

.700  x  10--' 

. 905  x  ID-2 

1  .02 
0. 9o 

Solut ion 

1  .1  39837 

(refined  nesh  and  Richardson  extrapolation) 

Case  2.  Increasing  Polvn 

orr.al  Order  on  2 A  Klenent  Me 

sh  ( No  Update  c f 

Hole 

Surface)  (Complete  Refinement) 

Element 

Type  (p) 

II -hit 

Error  Estimate 

ii e  Ik 

e 

Error  Estimate 

li e  1 E 

e* 

Linear  (1) 

l .139816 

.583  x  10_2 

.85 

.700  x  10-2 

1  .06 

Quadratic 

(2) 

1.159832 

.269  x  10-2 

.93 

.297  x  10"2 

1  .03 

Cubic  (3) 

1  . 1  39839 

.162  x  10-2 

•  9u 

.173  x  10-2 

1  .00 

"Exact " 

Solut ion 

1.1  59835  (refined  t.iesh  and  Richardson 

extrapol at  ion) 
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The  results  given  as  Case  2  in  the  table  are  for  a  uniform 
increase  of  the  polynomial  order  on  all  elements  based  on  the 
coarse  mesh  indicated  in  Figure  9.  Here  the  surface  geometry  of 
the  hole  has  not  been  updated  with  refinement.  Again,  the  error 
estimates  ||e||g  are  excellent  but  the  example  indicates  the 
practical  weakness  of  a  global  error  measure.  The  first  column 
in  the  table  indicates  that  the  significant  error  in  stresses  in 
the  immediate  region  of  the  hole  surface  appears  as  only  a  small 
perturbation  to  the  global  energy  of  the  region. 

Obviously  a  local  measure  is  required.  In  Figure  10  we  plot 
the  stress  tangential  to  the  hole  surface  in  two  elements  adja¬ 
cent  to  the  hole.  The  finite  element  stress,  plotted  as  a  linear 
interpolant  through  the  2x  2  Gauss  point  values,  is  scaled  using 
the  factors  FS^  given  by  (3.53),  The  order  of  the  error  in  the 
stresses  is  accurately  indicated  even  at  this  local  element  level 


1  EXACT  SOLUTION 
(INFINITE  REGION) 

2  FINITE  ELEMENT 
SOLUTION 

3  SCALED  BOUNDS 


Fig. 10.  Local  stress  errors:  _i  =  8b. -  on  2--elemont  mesh. 


4.  ADAPTIVITY 

The  concept  of  adaptivity  follows  naturally  from  the  previous 
discussion  as  the  expansion  of  the  trial  function  space  , 
hierarchically  or  non-hierarch ical Iv ,  but  only  where  the  space 
is  shown  to.  be  deficient.  It  has  been  taken  as  a  basis  for  adap¬ 
tivity  that  the  sequence  of  finite  element  solutions  must  follow 
the  best  rate  of  convergence  in  terms  of  the  number  of  degrees  of 
freedom  of  the  structure.  The  optimal  rates  of  convergence  for 
both  the  h  and  p  versions  of  the  finite  element  method  have 
been  identified  and  quantified  in  [13]. 

Two  programs  are  being  developed  for  this  research.  The 
first  is  based  on  the  h  convergence  process  in  Babuska  [7], 

as  the  indicator  for  mesh  refinement. 
i> 


which  utilizes  ||e||j.^ 
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The  second  p-convergence  algorithm  uses  directly  the  hierarchical 
elements  and  indicators  (3.38). 

The  strategy  for  selection  of  new  degrees  of  freedom  in 
adaptive  processes  is  not  uniquely  resolved.  Here  we  refine  on 
the  basis  of  evaluating  all  error  indicators  and  including  degrees 
of  freedom  whose  indicator  exceeds  one-half  of  the  maximum  value 
or,  in  the  case  of  the  h-convergence  program,  subdividing  elements 
whose  indicator  exceeds  this  value.  It  has  been  found  that  the 
path  followed  by  the  adaptive  process  is  not  greatly  affected  by 
changing  the  one-half  factor.  However,  in  practice  this  choice 
may  affect  the  expense  of  the  solution  process  and  alternative 
strategies  are  discussed  both  by  Babuska  [7]  and  Peano  [38], 

4.1  Examples 

4.1.1.  Cantilever  with  applied  end  moment.  The  error  estimators 
1 1  e  |  Ie  (£2  ^  are  identical  on  all  elements  in  the  first  two  meshes 

shown  in  Figure  8,  so  the  sequence  of  results  in  Table  1  corres¬ 
ponds  to  a  h-adaptive  process.  The  accuracy  of  the  error  esti¬ 
mators  indicates  chat  an  effective  stopping  criterion  on  the 
basis  of  the  energy  norm,  stress  or  displacement,  could  be  defined. 

4.1.2.  Small  hole  in  large  region.  Both  the  h  and  p  adaptive 
processes  have  been  applied  to  the  problem  defined  in  Figure  9 
and  the  results  plotted  on  Figure  11.  Notice  that  the  convergence 
is  expressed  in  terms  of  number  of  degrees  of  freedom  and  not  as 
usual  in  terms  of  h.  In  the  p-version  h  remains  constant  and 

a  comparison  would  not  be  possible.  For  2-D  elements  0(f>2)  = 

0(1/N)  so  we  expect,  for  example,  for  linear  elements,  slopes  £1. 

The  plots  in  Figure  11  can  be  divided  in  two  groups:  conver¬ 
gence  using  linear  trial  spaces  (adaptive  or  non-adaptive) ,  and 
convergence  using  higher  order  trial  spaces.  The  first  group 
includes  solution  extensions  nos.  3,  12,13,  14,  and  the  second 
group  extensions  nos.  2,4,  11,  from  meshes  Ml  and  M3.  Extensions 
4,  13,14  are  adaptive  and  the  type  of  meshes  obtained  are  repre¬ 
sented  in  Figures  12  and  13. 

As  expected,  the  rates  of  convergence  are  higher  for  the 
second  group  because  of  the  better  convergence  characteristics  of 
higher  order  finite  elements.  Within  each  group  the  adapaiive 
solutions  are  better  because  there  is  an  intelligent  criterion 
to  select  the  new  degrees  of  freedom. 

The  adaptive  p-extension  (extension  no. 4)  tends  in  the  limit 
to  the  cubic  solution  of  extension  no . 2 ,  due  to  the  fact  that  we 
set  a  limit  of  complete  cubic  modes  in  the  adaptive  p-convergence 
program.  This  last  result  shows  that  a  very  accurate  solution  on 
p-convergence  extensions  requi._~  the  use  of  higher  order  poly¬ 
nomials.  The  disadvantage  is  that  this  leads  to  a  loss  in  point- 
wise  convergence  because  of  the  'noise'  associated  with  the  poly¬ 
nomial  oscillations. 


A  mixed  h  and  p  convergence  model  could  therefore  be  the 
best  strategy.  In  [ 1A]  it  is  proved  that  the  combined  versions 
produce  higher  order  rates  of  convergence  than  either  the  h  or 
p-version  by  itself,  indicating  that  research  in  this  direction 
is  necessary.  We  note,  however,  that  the  resulting  program  will 
have  a  very  complex  structure.  In  this  context  we  reference  the 
work  done  on  the  h-version  [A3, 57]  which  by  itself  presents  a 
highly  complex  situation. 

The  p-convergence  programs,  on  the  other  hand,  have  a  simpler 
structure  but  nevertheless  more  elaborate  than  the  usual  finite 
element  codes,  since  we  have  to  allow  for  error  subroutines, 
automatic  node  generation  subroutines,  and  multi-level  finite 
element  types. 

The  meshes  obtained  for  the  h  and  p  adaptive  processes  of 
Extensions  1A  and  A  are  given  in  Figure  12  and  13  respectively. 
Similar  refinement  near  the  hole  is  seen  in  both  cases.  Finally, 
we  note  from  Table  3  that  the  error  estimators  for  the  h  adap¬ 
tive  process  given  in  Figure  12  are  again  accurate  enough  to  pro¬ 
vide  a  stopping  criterion.  However,  in  Section  A. 1.2  it  was 
seen  that  the  global  energy  norm  gave  little  indication  of  the 
accuracy  of  stresses  near  the  hole  surface.  The  local  stress 
error  estimator  discussed  in  that  section  may  provide  a  more 
practical  accuracy  test. 

Obviously  very  powerful  solution  algorithms  can  be  based  on 
these  processes.  The  efficiency  of  the  hierarchic  indicators  is 
best  shown  by  returning  to  the  cantilever  beam  of  Figure  8  and 
considering  a  tip  shear  load.  The  interpolants  required  on  the 
interfaces  between  elements  such  as  AB ,  in  the  figure,  depends 
on  the  Poisson's  ratio.  With  V  =0.0  only  cubic  interpolants 
in  the  x-direction  are  required;  with  V  =  0.3  the  exact  solution 
only  follows  if  quadratic  as  well  as  cubic  interpolants  are  added 
on  the  interfaces.  The  adaptive  process  based  on  the  hierarchic 
error  indicators  is  sensitive  to  exactly  these  requirements  and 
the  quadratics  are  left  out  of  the  adaptive  process  for  the  first 
problem. 
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5.  CONCLUDING  REMARKS 

A  general  description  of  'hierarchical'  finite  elements, 
including  both  h  and  p  versions  of  the  finite  element  method 
has  been  presented.  There  are  three  main  advantages  of  this 
approach.  Firstly,  it  leads  to  improved  conditioning  of  the 
stiffness  matrix  and  the  topology  of  the  stiffness  matrix  indi¬ 
cates  an  efficient  partitioning  for  block  iteration  solution 
procedures.  Secondly,  compatibility  is  easily  enforced  in 
meshes  with  a  graded  refinement  of  polynomial  order  or  element 
size.  Finally,  in  the  new  generation  of  adaptive  finite  element 
schemes  higher  order  hierarchical  modes  provide  an  indicator  for 
the  selection  of  the  new  degrees  of  freedom  which  should  be 
added  to  the  finite  element  mesh. 

We  have  also  demonstrated  in  this  paper  that  accurate  error 
estimators  are  available  and  can  be  interpreted  locally  at  least 
for  problems  with  stress  concentrations  no  greater  than  the  cir¬ 
cular  hole.  Projections  of  the  residuals  required  for  the  eval¬ 
uation  of  these  error  estimators  give  the  hierarchical  indicators 
which  control  the  adaptive  processes.  Naturally,  the  program 
architecture  becomes  complex  but  processes  allowing  the  accuracy 
required  of  the  solution  to  be  pre-specif ied  are  within  reach. 


APPENDIX 

PROOF  OF  THE  BOUND  OF  EQUATION  (3.30) 

The  energy  norm  of  the  error  is  given  by 
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